% Calculates GInd
%
% GInd is a vector of indexes. It describes how DIVIDING by the inflation 
% rate maps states to new states.
% 
% For example, if the value of the ith element of GInd is j, this means 
% that DIVIDING by inflation moves a firm in the ith state to the jth state
%
% A firm that is in a state (p_t-1(z)/P_t-1,A_t(z),S_t/P_t-1) will move to 
% state (p_t-1(z)/P_t,A_t(z),S_t/P_t) due to inflation between periods t-1 
% and t. If the initial state is the ith state. Then the ith element of 
% GInd gives the value of the state the firm goes to due to inflation.
%
% Jon Steinsson and Emi Nakamura, July 2006
%**********************************************************************

function GInd = CalculateGInd(a,rp,rad,G)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

ErVrpValue = repmat(rp,[1 agridsize radgridsize]);
ErVrpValue = permute(ErVrpValue,[2 1 3]);
ErVrpValue = reshape(ErVrpValue,[Ssize 1]);

ErVradValue = repmat(rad,[1 agridsize rpgridsize]);
ErVradValue = permute(ErVradValue,[2 3 1]);
ErVradValue = reshape(ErVradValue,[Ssize 1]);

G = repmat(G, [1 agridsize rpgridsize]);
G = permute(G,[2 3 1]);
G = reshape(G,[Ssize 1]);

ErVrpValue = ErVrpValue - G;
ErVradValue = ErVradValue - G;

ErVrpValue = (ErVrpValue < min(rp)*ones(Ssize,1)).*(min(rp)*ones(Ssize,1)) ...
    + (ErVrpValue > max(rp)*ones(Ssize,1)).*(max(rp)*ones(Ssize,1)) ...
    + (ErVrpValue >= min(rp)*ones(Ssize,1) & ErVrpValue <= max(rp)*ones(Ssize,1)).*ErVrpValue;

ErVradValue = (ErVradValue < min(rad)*ones(Ssize,1)).*(min(rad)*ones(Ssize,1)) ...
    + (ErVradValue > max(rad)*ones(Ssize,1)).*(max(rad)*ones(Ssize,1)) ...
    + (ErVradValue >= min(rad)*ones(Ssize,1) & ErVradValue <= max(rad)*ones(Ssize,1)).*ErVradValue;

ErVaInd = repmat((1:agridsize)',rpgridsize*radgridsize,1);

[dummy,ErVrpInd] = ismember(int32(ErVrpValue*10^6),int32(rp*10^6));
[dummy,ErVradInd] = ismember(int32(ErVradValue*10^6),int32(rad*10^6));
clear dummy
%ErVrpInd = match(int32(ErVrpValue*10^6),int32(rp*10^6));
%ErVradInd = match(int32(ErVradValue*10^6),int32(rad*10^6));
clear ErVrpValue ErVradValue

GInd = sub2ind([agridsize rpgridsize radgridsize],ErVaInd,ErVrpInd,ErVradInd);
